Hyperboloidal evolution of test fields in three spatial dimensions 
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We present the numerical implementation of a clean solution to the outer boundary and radiation 
extraction problems within the 3+1 formalism for hyperbolic partial differential equations on a given 
background. Our approach is based on compactification at null infinity in hyperboloidal scri fixing 
coordinates. We report numerical tests for the particular example of a scalar wave equation on 
Minkowski and Schwarzschild backgrounds. We address issues related to the implementation of the 
hyperboloidal approach for the Einstein equations, such as nonlinear source functions, matching, 
and evaluation of formally singular terms at null infinity. 



I. INTRODUCTION 

The most common method in the numerical construc- 
tion of solutions to hyperbolic partial differential equa- 
tions (PDEs) is the 3+1 formalism. In this approach the 
spacetime is foliated by a family of spacclikc hypcrsur- 
faces that are level sets of a time function. Typically, 
the spatially unbounded domain is truncated by an arti- 
ficial timelike outer boundary. The solution is then cal- 
culated numerically on a finite, spatially compact domain 
with given initial and boundary data. In this context, we 
present a successful numerical implementation of a clean 
solution to difficulties related to (i) the construction and 
imposition of outer boundary conditions, and (ii) the ex- 
traction of outgoing radiation for test fields. 

The first difficulty arises due to the artificial nature 
of the timelike boundary, which is not part of the phys- 
ical problem. Ideally, boundary conditions ensure trans- 
parency of the domain boundary; however, spurious re- 
flections occur from the outer boundary even for the sim- 
ple case of the flat wave equation in a three-dimensional 
ball pQ. Boundary conditions that are constructed to 
minimize spurious reflections are called nonreflecting, ab- 
sorbing or radiation boundary conditions [2H1]. They 
involve approximations to ideal conditions and substan- 
tial work goes into error controlling. In particular, it is 
difficult to construct, improve, or implement accurate ab- 
sorbing boundary conditions when the background curva- 
ture does not vanish or when non-linear terms appear in 
the equations. One needs to account for backscatter off of 
curvature or self-interaction of the field near the bound- 
ary. This is relevant because a bad choice of boundary 
data can reduce the accuracy and even destroy certain 
features of the solution. Especially in long time simula- 
tions, the outer boundary can be a significant limiting 
factor for the accuracy of numerical computations. 

The second difficulty is related to the notion of radia- 
tion. There are three scales of interest in radiative sys- 
tems: the scale of the source, the wavelength of the emit- 
ted radiation, and the location of the observer. Radiation 
is defined asymptotically in the far field zone where the 
observer is located. To study radiation accurately, one 
needs to evolve large portions of spacetime while keeping 
high resolution near the source. This is inefficient with 



Cauchy foliations because the resolution in the outer do- 
main is unnecessarily high even with mesh refinement 
techniques, and the time lag between the source and the 
observer increases proportionally to their distance. 

A clean solution to these difficulties is to avoid the ar- 
tificial truncation of the computational domain by com- 
pactification. It is well known that compactification at 
spatial infinity is not compatible with hyperbolic PDEs 
[5] . It leads to a loss of resolution due to blueshift in the 
wavenumber of outgoing waves. Compactification needs 
to be performed with respect to the outgoing character- 
istic direction of the hyperbolic PDE. In the 3+1 formal- 
ism a suitable time transformation needs to be applied to 
standard Cauchy surfaces before compactification so that 
the time surfaces approach null infinity in the asymptotic 
domain. Such time surfaces are called hyperboloidal be- 
cause their causal behavior resembles that of standard 
hyperboloids in Minkowski spacetime [BJ. 

Until now, numerical studies of the hyperboloidal ini- 
tial value problem for test fields have been restricted to 
one or two dimensional codes [7HT4] with the exception 
of [T5]. The hyperboloidal approach has been shown to 
be accurate and efficient in lower dimensional studies, 
but in |15j the authors report large numerical errors in 
the simulation of Maxwell fields across future null infin- 
ity in flat spacetime with a 3D code based on Cartesian 
coordinates. They deal with this problem by adding an 
artificial cosmological term to the background and by 
artificially tilting the light cones beyond null infinity to 
make grid boundaries spacelike. They suggest, however, 
that the errors may not be present in a numerical code 
that supports a spherical grid topology. We confirm this 
expectation in scri fixing coordinates with the Spectral 
Einstein Code SpEC |Tfj] . 

A strong motivation for the hyperboloidal approach 
is to improve the accuracy and efficiency of numerical 
relativistic calculations of binary black holes where the 
outer boundary and the radiation extraction problems 
have distinctive features [17-24 . There is an ongoing 
effort to solve a hyperboloidal initial value problem for 
various versions of the Einstein equations |25H31j . Our 
experiments have been chosen to address issues that ap- 
pear in these calculations. 
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II. HYPERBOLOIDAL SCRI FIXING 
COMPACTIFICATION 

It is favorable for numerical purposes to map null in- 
finity (scri) to a fixed coordinate sphere [3H 133] • Hy- 
perboloidal scri fixing coordinates for asymptotically flat 
spacetimes have been presented in [34] . Starting from 
the representation of a background spacetime based on 
Cauchy surfaces, hyperboloidal scri fixing compactifica- 
tion consists of: (i) a time transformation, (ii) a radial 
compactification, and (iii) a rescaling. 

The first class of hyperboloidal surfaces studied are 
constant mean curvature (CMC) surfaces [35^,, EE] • Brill, 
Cavallo and Isenberg constructed spherically symmet- 
ric CMC surfaces with nonvanishing mean extrinsic ex- 
trinsic curvature in Schwarzschild spacetime 37J. Re- 
cently, Malec and O'Murchadha analyzed such surfaces 
in great detail [3H1 [35] . CMC slicings serve as the start- 
ing point for the analysis of hyperboloidal slicings of 
the Schwarzschild spacetime defined by members of the 
Bona-Masso family of slicing conditions [40 . Beyond 
Schwarzschild spacetime, the CMC condition plays a cen- 
tral role in the construction of hyperboloidal initial data 
for Einstein equations presented by Buchman, Bardeen 
and Pfeiffer gT]. 

We employ CMC foliations in most of the numerical 
experiments. However, the CMC condition may be re- 
garded as restrictive due to its local nature, whereas 
the hyperboloidal condition for spacelike surfaces is only 
asymptotic. Therefore we discuss the steps in the con- 
struction of hyperboloidal scri fixing coordinates in a gen- 
eral manner for a given static background before present- 
ing the explicit results. 



A. General construction on a static background 

The metric on a static background spacetime can be 
written as 



We choose the function f2 such that it vanishes at a finite 
coordinate value with respect to p (see for example ( 10 ) 
or (14|). This implies that f2 vanishes at infinity with 
respect to r. The zero set of f2 corresponds to future null 
infinity denoted by (scri). 

The above compactification leads to a singular met- 
ric at the finite coordinate distance that corresponds to 
scri. This singularity can be rescaled away with a suit- 
able conformal factor if the spacetime is asymptotically 
simple |42j . We introduce the conformal metric 



9 = ti 2 g. 



(4) 



The metric g is not a solution to the vacuum Einstein 
equations due to the transformation behavior of the con- 
nection under a conformal rescaling of the metric [43j . 

The transformations Q, and Q lead to the con- 
formal metric 



g = Q 2 gtt dr 2 +2g tt H L drdp- 



~9ttH 2 



n 2 



L 2 dp 2 +p 2 da 2 . 

(5) 

We defined H := (dh/dr)(p) and L := 57 — pVl' . A prime 
denotes derivative with respect to p. Not every choice of 
the height function derivative H will result in a spacelike 
foliation. The height function derivative must satisfy 



H 2 < 



9r 



-gu 



(6) 



The formal singularity of the conformal metric ^ can be 
removed by asymptotic simplicity and by suitable choices 
of the involved functions as demonstrated explicitly for 
Minkowski and Schwarzschild spacetimes in the following 
subsections. 

We will often refer to the 3+1 decomposition of the 
conformal metric with respect to the time direction. In- 
troducing the lapse a, the shift /3, and the spatial metric 
function 7, we write 



g = (-a 2 + 7 2 /3 2 ) dr 2 + 2 7 2 /? dr dp + 7 2 dp 2 + p 2 da 2 



= gudt 2 + g Tr dr 2 + r 2 da 2 , 



(1) 



(7) 



where da 2 is the standard metric on the unit sphere. We 
introduce a new time function 



B. Minkowski spacetime 



t = t- h(r). 



(2) 



Here, h(r) is called the height function. This transfor- 
mation implies that the coordinates of the foliation are 
adapted to time symmetry in the sense that the time- 
like Killing vector is left invariant. We choose the height 
function so that level sets of r asymptotically approach 
future null infinity instead of spatial infinity (see for ex- 
ample ([9[> and pi])). 

The outgoing spatial direction is compactified by the 
radial transformation 



(3) 



fi(p)' 



The Minkowski metric fj in the standard spherical co- 
ordinate system reads 



fj = -dt 2 + dr 2 + r 2 da 2 . 



(8) 



We set the height function such that level sets of the new 
time coordinate r are CMC surfaces: 



/'('• = \/^ + r 2 , H(r)= J 

§2 + r 2 



(9) 



Here, a and S are free parameters. The mean extrinsic 
curvature of level sets of r reads K = 3S/a. Surfaces 
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with positive K approach future null infinity in our con- 
vention. The radial transformation is determined by the 
choice of the conformal factor. We set 



S 2 -p 2 
2a 



s 2 + P 2 

2a 



L 



(10) 



The set {p = S} corresponds to future null infinity and 
{/? = 0} corresponds to the origin. The conformal factor 
is a function of p 2 in Minkowski spacetime because the 
coordinate p is not regular at the origin. From ^ we get 
for the conformal Minkowski metric 



■q = n z fi = -fr dr 2 - 2p- dr dp + dp z + p z da 2 . (11) 

a 

Lapse, shift and the spatial metric function for the above 
metric are 



S 2 + P 2 R S 

a = — , p = -p— , 



2a ' a 

The Ricci scalar of the conformal metric reads 
12(S 2 - P 2 ){p 2 + 3S 2 ) 



7=1. (12) 



R 



(p 2 + S 2 ) 3 



C. Schwarzschild spacetime 

The Schwarzschild metric g with mass m in the stan- 
dard spherical coordinate system reads 

g = -Fdt 2 + —dr 2 + r 2 da 2 , with F = 1 - — . 

F r 



1. CMC foliation 

The height function for a CMC foliation of 
Schwarzschild spacetime can not be written in explicit 
form. Its derivative is given by |37H39j 



H = 



J 



fVJ 2 + f' 



J := 



Kr 



C 



(13) 



Here, K is the mean extrinsic curvature and C is an in- 
tegration constant. For numerical applications these pa- 
rameters should satisfy K > and C > 8m 3 K/3, which 
makes sure that the surfaces approach future null infinity 
and cross the future event horizon. 

A convenient choice for the conformal factor is 



n = l- 



s 



=> L = 1. 



(14) 



With J := flJ, we get for the conformal Schwarzschild 
metric 



g = -fl 2 Fdr 2 - 2 



J 



VJ 2 + Fn- 



■ dr dp 



1 

J 2 + Ffl 2 



dp 2 + p 2 da 2 . 



(15) 



The functions from the 3+1 decomposition ^ read 

a = \jj~ 2 + Ffi 2 , P = -Ja 7 = aT 1 . 

The Ricci scalar reads 

12fi(m(2p- S) + P S) 



(16) 



R 



p 2 S 3 



2. Matching 

The success of numerical calculations in computing 
the inspiral of binary black hole systems suggests that 
methods developed for the treatment of the asymptotic 
domain should be independent of the interior evolution. 
Also, many of the sophisticated techniques dealing with 
discontinuities and shocks related to the presence of mat- 
ter near the black holes use special techniques in specific 
coordinate systems. It is desirable to leave these meth- 
ods untouched in an effort to solve problems related to 
the asymptotic domain. 

An essential advantage of hyperboloidal foliations in 
this context is their flexibility. The only condition for a 
hyperboloidal foliation in compact domains is that it is 
spacelike. One can keep the numerical calculation near 
the sources untouched and change the foliation in the 
exterior domain to approach null infinity. In this set- 
ting, hyperboloidal coordinates are used only as a means 
to avoid boundary conditions and to have access to the 
idealized radiation signal [12l Q31 [34] . 

We describe the matching of an ingoing Eddington- 
Finkelstein foliation near a Schwarzschild black hole to 
hyperboloidal coordinates in an exterior domain. The 
Schwarzschild spacetime with an arbitrary H and can 
be written as 



g = -fl 2 Fdr 2 -2L 
L 2 ( „ 2mQ, 



2mD, 
P 



(1 + H)) drdp 



1 - H z 



(l + H) 2 ) dp z +p z da z . (17) 



n 2 \ p 

The functions of the 3+1 decomposition are 



l-H 2 



2m!l 



(l + H) 2 



7 = 



(18) 



We recover the Schwarzschild metric in ingoing 
Eddington-Finkelstein coordinates for f2j n = 1 and H ln = 
0. We use this choice in an inner domain p < p ln . To get 
hyperboloidal scri fixing coordinates outside a compact 
domain we set onp> p out (see |34j ) 



On 



1 



S' 



1 



4mfi (8m 2 - C^)fl 2 



S 



(19) 
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The matching between the domains is performed by a 
smooth transition function 



0, P < Pin, 
f = { fa, Pin < P < Pout, 

1, P > Pout, 



(20) 



where f T is O [45 



11 / s 
fr = — I — tanh — ( tan w(p) 
2 2 V 7T 1 



tantti(p) 



with 



w(p) :-- 



7T p - P h 



2 Pout — Ph 



The free parameter q determines the point p\i% at which 
fx takes the value 0.5 and the parameter s determines 
the slope of fr at this point. 



We set the functions in ( 18 1 as follows 



n = i-/£ L = n- P n', H = fH b 



The Ricci scalar for the metric (17) reads 
6fl 



R 



L 3 r 2 



{2LQ!{m9. + mrQ! - r) - VLr(r - 2mf2)f2")- 



III. THE CONFORMAL METHOD 

We take the wave equation with a source term as a toy 
problem 



(21) 



The scalar field ^> is evolved in a background spacetime 
with metric g. Here, Dg :— g^V^V,, where V is the 
Levi-Civita connection of the metric g. 

The scalar wave equation is invariant under a change 
of coordinates but not under a conformal rescaling of 
the background metric. We need to take the conformal 
transformation behavior of the scalar wave equation into 
account to solve it on a conformal background. A con- 
formal rescaling g = fl 2 g of the physical metric g with a 
conformal factor Q > implies the transformation 



□ - -r) v = q- 3 

6 



□ - 1 *. 



with * = 



where R and R are the Ricci scalars of the rescaled and 
physical metrics g and g respectively. The scalar wave 
equation ( 21 ) on the background of a solution of the vac- 



uum Einstein equations is equivalent to 



1 



- -R^ = n-^F(x' i , Q-'g, 



(22) 



This equation is singular only if the source term does 
not fall off sufficiently fast towards infinity. This is not a 



strong restriction because physically motivated problems 
typically have source terms that have strong fall off or 
compact support. 

The conformal transformation is a minor modification 
that does not change the principal part of the system. It 
introduces only an additional source term to the scalar 
wave equation. The results we present in this paper are 
expected to apply similarly for any system of equations 
for test fields with a reasonable conformal transformation 
behavior. For example, Maxwell and Yang-Mills fields 
are conformally invariant and therefore, the application 
of the hyperboloidal method does not require modifica- 
tions of the geometric equations. The advantage of such 
a general geometric approach is that, in contrast to the 
construction of artificial outer boundary conditions, a dif- 
ferent numerical technique, background, or equation does 
not significantly change the application of the method. 
The situation is rather different for the Einstein equa- 
tions where the metric plays a double role of background 
and unknown variable 1461. 



IV. NUMERICAL EXPERIMENTS 



A. Implementation 

We solve the conformal wave equation with the Spec- 
tral Einstein Code SpEC [TB] used for binary black hole 
evolutions (13 HI]- A pseudospectral technique is applied 
on each subdomain to evolve the scalar wave equation in 
time. An important property of the numerical calcula- 
tion is that it employs spherical boundaries. The numer- 
ical outer boundary of the computational domain corre- 
sponds to future null infinity, whose cuts are topological 
spheres, therefore it is advantageous to use spherical grid 
topology. Our basis functions are Chebyshev polynomi- 
als in the radial direction and spherical harmonics in the 
angular directions. 

Numerical codes that work with multidomain tech- 
niques typically employ a global Cartesian system. We 
give for convenience the transformation of a spherically 
symmetric metric of the form ([7| into Cartesian coordi- 
nates. We have coordinates Xi such that p 2 — SijXiXj, 
where is the 3D Kronecker delta. The 3+1 decom- 
position of a metric g in Cartesian coordinates can be 
written as 



g = (-a 2 + ^P l P 3 ) dr 2 + 2^0* dr dx j + 7ij dx l dx j . 



We have 



P P 
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The derivatives are calculated by 

2 



X^Xj iCfa 



p° 



T 



1 



p \ p 

(3\ XiX k 



P 



P 



ah 



We solve a hyperboloidal initial value problem for the 



conformal wave equation (22 1 in first order form with 
auxiliary variables $i := ai^> and IT := — l/a(d T ^ — 
p l di^) as in ED] ■ As initial data we set all fields to 
zero except 11(0, p) which is set to an off-centered Gaus- 
sian pulse 



n(o,p) = e -i x - x °i 1° , 

or to a pure Yj m mode 

n(o, p) = e-('-'°> a / ff2 y, m (i? )¥ o. 

B. Minkowski spacetime 

1. The setup 



(23) 



(24) 



problem. It simplifies the implementation and makes er- 
ror controlling and stability considerations trivial. 

Our simulation domain in the radial direction is given 
by p £ [0,20]. The domain structure includes a cube 
around the origin with the domain Xi £ [—2,2] and 4 
spherical shells extending to future null infinity (Fig. [TJ . 




Figure 1: The numerical grid topology for calculations in 
Minkowski spacetime. We have a cube around the origin with 
the domain [—2, 2] in each Cartesian direction and 4 spherical 
shells extending from p — 2 to future null infinity at p = 20. 
The colors depict an off-centered Gaussian initial data for II 
as given in ( 23 1 . 



The gauge parameters in our coordinatization of the 
conformal Minkowski spacetime (111 are the coordinate 



location of null infinity S, and a parameter a that con- 
trols the mean extrinsic curvature of the hyperboloidal 
surfaces K via K = 3S/a. To set these parameters suit- 
ably we discuss the characteristic structure of Minkowski 
spacetime in our gauge. The coordinate speeds of in- and 
outgoing radial light rays are given by 



c±=±~ 



We get with (12 1 



i 

2a 



(p-S? 



(25) 



(26) 



The outgoing speed near y + is proportional to S 2 /a. 
We need to keep the characteristic speeds at the order 
of unity to avoid strong restrictions on the time step 
due to the Courant-Friedrichs-Lewy (CFL) condition. 
We choose the parameters such that a ~ S 2 implying 
K ~ 1/S. In our calculations we set S — 20 and a = 400. 



The formulae ( 26 ) show that there are no incoming 
characteristics from the outer boundary at p — S imply- 
ing that no boundary conditions are needed or allowed 
at the numerical outer boundary. This property of hy- 
perboloidal scri fixing compactification is the main ad- 
vantage over standard treatments of the outer boundary 



2. Convergence 

Spatial truncation errors converge exponentially in a 
pscudospectral code. We show such spectral convergence 
in Fig. [2] for an evolution with vanishing source and off- 
centered initial data (Eq. ( 23 ) and Fig. [l]) . We plot the 
L2-norm of the constraint field 



(27) 



in time. The convergence properties of pure Yi m -data 
are similar. The off-centered case provides a stringent 
test on the accuracy of the code because the evolution is 
highly asymmetrical and the resolution requirements are 
higher. 



3. Semilinear wave equation 

The source-free wave equation on Minkowski space- 
time is a simple test because the wave package leaves 
the spacetime along the characteristic directions due to 
the Huygen's principle. When nonlinear source terms are 
present, however, there is backscatter due to the self in- 
teraction of the field and the numerical method needs 
to deal with nonlinearities. This is a strong test for 
transparency boundary conditions because they must not 
eliminate all reflections from the outer boundary but only 
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Figure 2: Spectral convergence for the evolution of off- 
centered initial data in Minkowski spacetime by the L2 norm 
of the constraint field d ( 27 1 . 



spurious ones because backscatter near the boundary is 
part of the solution. 

We demonstrate the applicability of our method to 
wave equations with nonlinear source terms by taking 
a source of the form F = — <F P . The minus sign corre- 
sponds to focusing 51 . The conformal wave equation 
with such a power source term reads 

- -m = -n p - 3 

6 

The right hand side is regular for p > 3. A good measure 
for the accuracy of the code is whether the backscatter 
of the field due to self-interaction can be calculated cor- 
rectly in long time evolutions. We show in Fig.[3]the local 
decay rates of the field for p — 3 as measured by various 
observers. The decay rates at the end of our simulation 
are —1.02 along null infinity and —1.97 along r = 30, 
very close to the expected asymptotic rates of —1 along 
null infinity and —2 along timelike surfaces. 

Decay Rates 
-1.0 




200 400 600 800 1000 



Figure 3: Local decay rates measured by various observers for 
the cubic wave equation. The locations of the observers vary 
from null infinity (top curve) to r = 30 (bottom curve). 

The implementation of nonlinear source terms poses 



no difficulties and does not change the numerical outer 
boundary treatment in the hyperboloidal method as long 
as the terms have sufficiently fast fall off. This behavior 
should be contrasted to accurate transparency boundary 
conditions that need to be modified depending on the 
details of the source terms. 



C. Schwarzschild spacetime 

1. The setup 

The gauge parameters in our coordinatization of the 
conformal Schwarzschild spacetime in a CMC foliation 
(15) are the coordinate location of null infinity, S 1 , the 



mean extrinsic curvature of the hyperboloidal surfaces, 
K, and an integration parameter C that influences the 
causal properties of the CMC foliation. A useful discus- 
sion of these parameters can be found in |41j . 



Evaluation of the characteristic speeds ( 25 1 at fu- 



ture null infinity shows that the outgoing speed becomes 
^K 2 S 2 . Therefore we choose the mean extrinsic curva- 
ture as K ~ 1/S 2 . We set the coordinate location of 
future null infinity to S = 20 as for Minkowski spacetime 
in the previous section. We further set m — 1, K = 0.07 
and (7=1. The event horizon is located at 



Ph 



2mS 
2m + S' 



We use the excision technique by placing the inner 
boundary inside the event horizon, where all characteris- 
tic fields are outgoing and no boundary conditions need 
to be applied. There is a coordinate singularity inside 
the event horizon for CMC surfaces on Schwarzschild 
spacetime given by the zero set of the lapse function [37l - 
|3"5I HTj . For the above choice of parameters the singular- 
ity is at p s = 1.732 and the horizon is at pu = 1.818. 
To avoid the singularity we put the numerical inner 
boundary just inside and very close to the event hori- 
zon. Our simulation domain in the radial direction reads 
p € [1.8,20]. The coordinate speeds of in- and outgoing 
radial characteristics on this domain with the above pa- 
rameters are plotted in Fig. ij The ingoing speed (red 
curve) vanishes at the outer boundary, whereas the out- 
going speed (blue curve) is non-vanishing and less than 
unity. The outgoing speed is fairly constant over the 
grid, which implies even resolution for the outgoing wave 
across the numerical domain that covers an infinite phys- 
ical domain. There are no boundary conditions needed 
because no characteristics enter the simulation domain. 
The domain structure consists of 9 spherical shells that 
extend from inside the event horizon at p = 1.8 to future 
null infinity at p = 20. 
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Figure 4: Coordinate speeds of ingoing (red) and outgoing 
(blue) radial characteristics in Schwarzschild spacetime with 
respect to a CMC foliation. There are no incoming charac- 
teristics to the simulation domain, so no boundary conditions 
need to be applied. 



2. Convergence 

We show spectral convergence in Fig. [5] for an evo- 
lution with vanishing source and off-centered data in 
Schwarzschild spacetime. For this plot we cover the sim- 
ulation domain with 9 spherical shells and fix the angular 
resolution. We perform the radial convergence test with 
N radial collocation points in each shell. Fig. [5] shows 
the L2 norm of the constraint field Cj defined in ( [27] ) 
over the whole grid including the outer boundary. 
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Figure 5: Convergence for the evolution of off-centered ini- 
tial data in Schwar zsch ild spacetime by the L2 norm of the 
constraint field d (27 1. The domain is covered by 9 spher- 
ical shells with radial N collocation points each. We fix the 
angular resolution. 



3. Quasmormal mode ringing and tail decay rates 

A perturbation of Schwarzschild spacetime evolves in 
three phases: an initial transient phase that depends on 
initial data, a quasinormal mode ringing, and a poly- 



nomial decay. We show the quasinormal mode ringing 
and the subsequent polynomial decay in Fig. [6] for the 
evolution of a pure I = 2, m = initial data in a half- 
logarithmic plot. The ringing has the form 



*(r) = a e"" 2T smfar + if). 



(28) 



Here, u>i and uj^ are the mode frequencies, a is the 
amplitude and ip is the phase of the wave signal. We 
measure these parameters by fitting the signal to the 
above formula using a simple least squares method on 
the interval r <E [60m, 120m]. We find bj x = 0.48389 
and U2 = 0.09661. These numerical values are very 
close to those obtained by using Leaver's continued frac- 
tion method [52j [53], which read u>i — 0.48364 and 
uj 2 = 0.09676. The fitting error in the mode frequen- 
cies dominate the numerical discretization error of the 
solution. 




100 200 300 400 500 

Figure 6: The quasinormal ringing and subsequent polyno- 
mial decay in Schwarzschild spacetime for initial data with 
Z = 2,m = 0ina half-logarithmic plot. 



In general, it is difficult to calculate the polynomial 
decay rates accurately with a 3D code. Furthermore, nu- 
merical calculations of such decay rates usually measure 
the rates near the black hole, which are different from the 
rates in the far-field zone [S3]. The correct and accurate 
numerical calculation of asymptotic decay rates as mea- 
sured by idealized observers can only be achieved in com- 
pactified schemes that include null infinity [TTJ [SH [S5J . 

We show in Fig. [7] the local power indices as measured 
by observers ranging from near the black hole to future 
null infinity for initial data with a pure I = 2, m = 2 
mode, which is expected to be the dominant mode for 
gravitational waves. This mode decays asymptotically in 
time with rates of —4 along null infinity and —7 along 
timelike surfaces. Fig. [7] indicates that these asymptotic 
decay rates can be reproduced cleanly in a 3D long time 
evolution with the hyperboloidal method. The coordi- 
nate location of the numerical outer boundary is 20. In 
the standard approach this would imply that Fig.[7]shows 
the decay rates for 50 crossing times without any contam- 
ination from the outer boundary. Of course, the notion 



of a crossing time does not apply in our case because 
the incoming propagation speed from the outer bound- 
ary vanishes. 

Decay Rates 



—4.0 c 




T 



Figure 7: Local decay rates measured by various observers for 
the evolution of pure / = 2, m — 2 initial data. The locations 
of the observers vary from null infinity (top curve) to r = 30m 
(bottom curve). 
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Figure 8: Convergence for the evolution of a pure I = 2, m = 2 
initial data in Schwarzschild spacetime with a matched folia- 
tion. 



4- Matching 

The matching is a rather arbitrary way of gluing hy- 
persurfaces together; therefore there are many free pa- 
rameters to be set. We fix the coordinate location of 
null infinity to S — 20 as in the previous computations. 
To set the constant Ct in ( 19 ) we consider the outgoing 
characteristic speed at null infinity, which reads S 2 /C T . 
So we set Ct = S. 



The formula for the height function derivative ( 19 ) can- 



not be used everywhere because it violates the spacelike 
condition With m = 1 and the above choice of pa- 
rameters, the condition is violated at r = 12.5. The in- 
ner matching radius must be larger than this value. We 
can move this point closer to the black hole in terms of 
the grid coordinate by performing the compactification 
earlier than the bending up of the height function. To 
achieve this we introduce parameters pf^, p^ ut , p?, p^ ut so 
that the matching of the physical coordinate to the com- 
pactifying coordinate is performed independently than 
the matching of the ingoing Eddington-Finkelstein time 
surfaces to asymptotically hyperboloidal time surfaces. 
A good choice of parameters that we used for the evolu- 
tions underlying the convergence plot in Fig. IS] is pP = 



1 Pout 



10, p H = 10, p H = 15. 



Fig. [8] shows the convergence properties of matched 
evolutions for initial data with a pure I = 2, m = 2 
mode. On the one hand, it is reassuring that match- 
ing can be performed stably in combination with pseu- 
dospectral methods. On the other hand, Fig. [8] shows 
that the errors are about four orders of magnitude big- 
ger than the errors in Fig. [5] so the matching is very 
inefficient compared to CMC calculations. This differ- 
ence can be explained as follows. The spatial redshift 
effect PJjl H31 GH] leads to highly efficient numerical cal- 



culations in combination with high order discretization 
schemes such as pseudospectral methods as we demon- 
strated in the previous sections. With matching, this 
effect is active only in part of the domain. Instead, the 
additional structure from matching needs to be resolved 
and is a source of numerical error. One can expect that 
lower order or finite element type numerical methods will 
work better with matching. 

The form of the metric that we implement for match- 
ing is formally singular at future null infinity where the 
conformal factor Q vanishes ( [18] ) . Instead of applying the 
analytic limits at the outer boundary, we employ in the 
outermost shell Gauss-Radau collocation points, which 
do not involve the point at the outer boundary. This fa- 
cilitates a natural extrapolation within the pseudospec- 
tral method and allows us to evolve the wave equation 
without numerically evaluating the formally singular co- 
efficients. A similar formal singularity is present in the 
conformal Einstein equations. A strategy to implement 
the hyperboloidal approach in combination with the gen- 
eralized harmonic formulation of the conformal Einstein 
equations requires the numerical evaluation of such for- 
mally singular terms at the outer boundary 29 . Our 
results suggest that a similar extrapolation method for 
the conformal Einstein equations may be successful. We 
emphasize, however, that the model problem studied here 
is simpler due to the decoupling of the formal singularity 
from the evolution equations. A successful implementa- 
tion of conformal Einstein equations that does not require 
the evaluation of such terms has been presented by Rhine 
using a constrained scheme in axisymmctry |31j . 



V. CONCLUSIONS 

The hyperboloidal method allows us to numerically 
calculate solutions to hyperbolic PDEs on unbounded do- 
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mains thereby avoiding the outer boundary problem and 
cleanly solving the radiation extraction problem. The 
approach in this paper is based on the conformal method 
[55] . the hyperboloidal initial value problem [3] and the 
hyperboloidal scri fixing compactification of asymptoti- 
cally flat spacetimes [34] in combination with the Spectral 
Einstein Code SpEC developed by the Caltech-Cornell 
collaboration [TB] . Our results obtained for the toy model 
of a scalar field can be expected to hold equally for other 
test fields such as Maxwell or Yang-Mills fields. Nu- 
merical tests with off centered and pure spherical har- 
monic data and nonlinear source terms in Minkowski and 
Schwarzschild spacetimes support our conclusions. 

The efficiency of the hyperboloidal method in combi- 
nation with pseudospectral methods is striking. All com- 
putations for this paper have been performed on a laptop 
without parallelization. Decay rates as strong as —7 can 
be accurately calculated in 3D with less than 200 ra- 
dial collocation points to cover an infinite domain of the 
Schwarzschild spacetime. 

There are many advantages of the hyperboloidal 
method beyond efficiency, especially in comparison to 
standard outer boundary treatments. First of all, the 
question of well-posedness and stability of boundary 
conditions becomes trivial because there are no bound- 
ary conditions. In standard approaches, boundary data 
needs to be chosen to ensure the transparency of the 
boundary and to model the exterior domain. This calcu- 
lation is typically complicated and falls out completely in 
the hyperboloidal method. The geometric nature of the 
method makes it largely independent of the details of nu- 
merical schemes, whereas it is a nontrivial task to imple- 
ment certain accurate boundary conditions for high order 
numerical discretizations. The hyperboloidal method re- 
quires minor modifications for a large class of equations 
and background spacetimes. Achieving good accuracy for 
the boundary treatment of each of the examples we pre- 
sented took decades of work and further improvements 
are exceedingly complicated, whereas the hyperboloidal 
method applies easily to problems with a suitable asymp- 
totic behavior. The method is also practical. There is no 
overhead in software implementation and therefore no 
need to consider runtime cost at the boundary in com- 
parison to the interior grid work. There is no need for 
error controlling because there are no errors other than 
numerical discretization. In addition, we have access to 
the radiation signal as measured by idealized observers 
at future null infinity. 

The time transformation can be regarded as a disad- 
vantage of the hyperboloidal method because certain ap- 
plications require a specific coordinate system in a com- 
pact domain. The matching technique avoids the time 
transformation in such a compact domain. We observe, 
however, that matching is less efficient than a CMC fo- 
liation. This is in accordance with expectations because 
with matching there is additional structure that needs 
to be resolved by the numerical scheme. It also intro- 
duces many free parameters complicating the application 
of the method. In cases where the use of a specific coordi- 



nate system is necessary, the decision whether to employ 
the matching technique will depend on the importance 
of boundary conditions, time scale of the simulation, and 
even the numerical scheme. Further research in this area 
may show that matching works well with certain numer- 
ical schemes, but for now our experiments suggest that 
matching should be avoided if possible. 

Currently, matching is required for the application of 
the hyperboloidal method in Kerr spacetime. A study 
of CMC foliations of Kerr spacetime would be beneficial 
for improving the accuracy and efficiency of numerical 
computations on a rotating black hole background. The 
analogous problem for nonrotating black holes has been 
solved in 1980 by Brill, Cavallo, and Isenberg [37]. It 
seems that time has come to attack the problem for ro- 
tating black holes. 

A further disadvantage of the hyperboloidal method is 
the requirement of a spherical grid topology at the outer 
boundary. Cuts of future null infinity naturally have 
spherical topology. Previous studies of hyperboloidal 
compactification with Cartesian grids report spurious re- 
flections of outgoing signals through null infinity |15j . 
Therefore this technical requirement may be a limiting 
factor for applying the hyperboloidal method in numeri- 
cal codes employing Cartesian grids. 

Eventually, the main problem that we would like to 
tackle is the hyperboloidal initial value problem in gen- 
eral relativity. Our results should be useful both for reg- 
ular conformal field equations 25 27J and for conformal 
Einstein equations [2"9Tf3"T] . Specifically relevant for the 
full problem in general relativity are the following ob- 
servations: free evolution is very accurate and efficient 
so that even strong decay rates can be resolved in 3D 
without using large computational resources; the spa- 
tial redshift effect of hyperboloidal foliations works well 
with pseudospectral methods; nonlinear source terms can 
be handled without difficulties; the coordinatization of a 
compact domain can be chosen freely but matching is an 
additional source of numerical error; and the numerical 
evaluation of formally singular terms at the outer bound- 
ary using Gauss-Radau collocation points leads to stable 
evolutions. These results may play a fundamental role in 
improving the accuracy and efficiency of numerical rela- 
tivistic calculations with the hyperboloidal approach. 
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